#include <iostream>
#include "spline.h"
#include <cmath>
#include <vector>
#include <algorithm>
#include <fstream>

class Fun1 : public Function//右上部分
{
public:
    double operator()(double _x)
    {
        return 2.0/3*pow(3-_x*_x, 0.5)+2.0/3*pow(_x, 0.5);
    }
    double diff(double _x)
    {
        return -2.0/3*_x*pow(3-_x*_x, -0.5)+1.0/3*pow(_x, -0.5);
    }
    double diff_2(double _x)
    {
        return -2.0/3*pow(3-_x*_x, -0.5)-2.0/3*_x*_x*pow(3-_x*_x, -1.5)-1.0/6*pow(_x, -1.5);
    }
};
class Fun2 : public Function//右下部分
{
public:
    double operator()(double _x)
    {
        return -2.0/3*pow(3-_x*_x, 0.5)+2.0/3*pow(_x, 0.5);
    }
    double diff(double _x)
    {
        return 2.0/3*_x*pow(3-_x*_x, -0.5)+1.0/3*pow(_x, -0.5);
    }
    double diff_2(double _x)
    {
        return 2.0/3*pow(3-_x*_x, -0.5)+2.0/3*_x*_x*pow(3-_x*_x, -1.5)-1.0/6*pow(_x, -1.5);
    }
};

class Discrete_function : public Function
{
private:
    std::vector<double> Point, value;
    int n;
public:
    Discrete_function(std::vector<double> _Point, std::vector<double> _value):Point(_Point), value(_value) 
    {
        n = Point.size();
    }

    double operator()(double _x)
    {
        for (int i = 0; i < n; i++)
        {
            if (Point[i] == _x)
                return value[i]; 
        }
        return 0;
    }

    double diff(double _x)
    {
        return 0;
    }

    double diff_2(double _x)
    {
        return 0;
    }
};

int main()
{
    std::ofstream fout1_1, fout1_2;
    fout1_1.open("E_point3_1.txt");
    fout1_2.open("E_point3_2.txt");
    double a = 0;
    double b = pow(3, 0.5);
    int N = 160;
    int n = 1000;
    int order = 3;
    int way = 3;
    Fun1 f1;
    Fun2 f2;
    double ccl;//累积弦长
    double length;
    std::vector<double> Point1;
    std::vector<double> Point2;
    std::vector<double> Point_x;
    std::vector<double> Point_y1;
    std::vector<double> Point_y2;
    ccl = 0;
    Point1.push_back(0.0);
    for(int j = 1; j <= N; j++)
    {
        ccl = ((b-a)/N)*((b-a)/N) + (f1(a + j*(b-a)/N) - f1(a + (j-1)*(b-a)/N))*(f1(a + j*(b-a)/N) - f1(a + (j-1)*(b-a)/N));
        ccl = sqrt(ccl);
        Point1.push_back(Point1[j-1]+ccl);
    }
    for(int j = 0; j <= N; j++)
    {
        Point_x.push_back(a + j*(b-a)/N);
        Point_y1.push_back(f1(a+j*(b-a)/N));
    }
    Discrete_function F1_x(Point1, Point_x);
    Discrete_function F1_y(Point1, Point_y1);
    B_spline_interpolation P1_x(F1_x, Point1, way, order);
    P1_x.solve();
    length = 0;
    while(length <= Point1[N])
    {
        fout1_1 << P1_x(length) << std::endl;
        length = length + 0.01;
    }
    B_spline_interpolation P1_y(F1_y, Point1, way, order);
    P1_y.solve();
    length = 0;
    while(length <= Point1[N])
    {
        fout1_1 << P1_y(length) << std::endl;
        length = length + 0.01;
    }

    Point2.push_back(0.0);
    for(int j = 1; j <= N; j++)
    {
        ccl = ((b-a)/N)*((b-a)/N) + (f2(a + j*(b-a)/N) - f2(a + (j-1)*(b-a)/N))*(f2(a + j*(b-a)/N) - f2(a + (j-1)*(b-a)/N));
        ccl = sqrt(ccl);
        Point2.push_back(Point2[j-1]+ccl);
    }
    for(int j = 0; j <= N; j++)
    {
        Point_y2.push_back(f2(a+j*(b-a)/N));
    }
    Discrete_function F2_y(Point2, Point_y2); 
    Discrete_function F2_x(Point2, Point_x); 
    B_spline_interpolation P2_x(F2_x, Point2, way, order);
    P2_x.solve();
    length = 0;
    while(length <= Point2[N])
    {
        fout1_2 << P2_x(length) << std::endl;
        length = length + 0.01;
    }
    B_spline_interpolation P2_y(F2_y, Point2, way, order);
    P2_y.solve();
    length = 0;
    while(length <= Point2[N])
    {
        fout1_2 << P2_y(length) << std::endl;
        length = length + 0.01;
    }
    fout1_1.close();
    fout1_2.close();
    return 0;
}